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Abstract: We compute the rapidity distributions of W and Z bosons produced at the 
Tevatron and the LHC through next-to-next-to leading order in QCD. Our results demon- 
strate remarkable stability with respect to variations of the factorization and renormaliza- 
tion scales for all values of rapidity accessible in current and future experiments. These 
processes are therefore "gold-plated": current theoretical knowledge yields QCD predic- 
tions accurate to better than one percent. These results strengthen the proposal to use W 
and Z production to determine parton-parton luminosities and constrain parton distribu- 
tion functions at the LHC. For example, LHC data should easily be able to distinguish the 
central parton distribution fit obtained by MRST from that obtained by Alekhin. 
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1. Introduction 



Drell-Yan production of lepton pairs through electroweak (EW) gauge bosons at hadron 
colhders occupies a special place in elementary particle physics. Historically, the Drell- 
Yan mechanism [1] was the first application of parton model ideas beyond deep inelastic 
scattering, and was later the route to discovery of the W and Z bosons [2]. Currently, it 
provides precise determinations of several Standard Model (SM) parameters, and places 
stringent constraints on many forms of new physics. Studies of W production at the 
Tevatron lead to determinations of the mass and width of the W boson with precision 
competitive with LEP2 measurements [3, 4, 5, 6]. The ratio of production cross sections 
for W and Z bosons, weighted by their leptonic branching fractions, is very accurately 
predicted in the Standard Model, and has been studied extensively at the Tevatron [7]. The 
rapidity distribution for produced Z bosons [8], and the charge asymmetry in leptons from 
W production [9], have also been measured at the Tevatron; both distributions are sensitive 
to the distribution of partons within the proton. Searches for non-standard contributions to 
the production rate of lepton pairs with invariant masses larger than Mw,z can be used to 
detect additional gauge bosons, such as the Z states that appear in most extensions of the 
SM. More generally, these searches constrain possible contact interactions between quarks 
and leptons arising from new physics at energy scales beyond those currently accessible 
[10]. 

With Run II of the Tevatron producing data, and with the LHC scheduled to begin 
operation shortly, an enormous number of W and Z bosons will soon be collected. This 
will significantly increase the precision of electroweak measurements, and will dramatically 
boost the sensitivity of new physics searches. To fully utilize these results, precise theoret- 
ical predictions for W and Z cross sections are needed. Current calculations are limited 
by uncertainties in parton distribution functions, as well as higher-order QCD and EW 
radiative corrections [11]. 

Parton distribution functions (PDFs) are determined from a global fit to a variety of 
data; unfortunately there is no direct experimental information for the combined values of 
(10^ GeV^) and Bjorken x (10"'^ to 10~^) that are relevant for electroweak physics at 
the LHC. PDFs for these parameter values are obtained through perturbative evolution 
of fits to PDFs at lower values of Q^, using the DGLAP equation. The complete results 
for the DGLAP evolution kernels at next-to-next-to leading order (NNLO) are not yet 
available. An approximate set of evolution kernels is used instead [12]. 

There are currently two sets of PDFs extracted with NNLO precision, using these 
approximate kernels. The MRST set [13] utilizes a broad variety of data; the drawback of 
this procedure is that the data set includes observables for which NNLO QCD corrections 
are not known. Alekhin's PDFs [14] are based on deep-inelastic scattering data only; 
this data set is somewhat restricted, but higher order QCD corrections can be included 
consistently [15]. The two PDF sets lead to slightly different (at the few percent level) 
predictions for the total rate of W and Z production at the Tevatron and the LHC [16, 
17]. We take this difference as a rough estimate of the current uncertainties in the PDFs 
needed for Tevatron and LHC physics; the individual PDF fits now also contain intrinsic 
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uncertainty estimates [14, 16]. 

The QCD corrections to EW gauge boson production have been studied by several 
groups. The complete NNLO corrections to the total cross section were computed some 
time ago [18, 19]. However, the total cross section is not an experimental observable, and 
significant extrapolations are required to compare this prediction to experiment. Ideally, 
one would like an event generator, at least at the parton level, which retains the full 
kinematics of the process and incorporates higher-order radiative corrections. Although 
there has been some recent progress towards this goal at NNLO [20, 21, 22, 23, 24], its 
completion will probably not occur for some time. 

NLO QCD corrections to more differential quantities in EW gauge boson production, 
including the vector boson rapidity distribution, were computed in Ref. [25]. A general- 
ization of this result to W and Z production at the Tevatron and the LHC yields NLO 
corrections of approximately 20% to 50%, and scale variations of a few percent. Since the 
NLO corrections are rather large, while the residual scale dependence is small, the actual 
reliability of the NLO results has been somewhat unclear. Even taking the NLO scale 
variations seriously, it is apparent that our knowledge of higher-order QCD corrections to 
EW gauge boson production is accurate to at best a few percent. 

This few percent precision in our knowledge of both PDFs and radiative corrections 
must be compared to the needs of the Tevatron and LHC physics programs. The W 
mass should be measured with a precision of ±30 MeV during Run II in each Tevatron 
experiment [26]; this uncertainty will be further decreased to ±15 MeV at the LHC [27]. 
Such a measurement strengthens the constraints that the precision EW data imposes on 
the Higgs mass and on many indirect manifestations of new physics. A precise theoretical 
prediction for M]^ requires knowledge of the W transverse momentum spectrum, as well 
as a good understanding of the relevant PDFs. To calibrate the detector response for the 
measurement of the W decay products, both the rapidity and p± spectra of the Z must 
also be well understood [26] 

At the LHC, many additional measurements will also require theoretical predictions 
accurate to a percent or better. The extremely large cross section for the Drell-Yan pro- 
cess at the LHC allows measurements of the Z boson rapidity distribution, and of the 
pseudorapidity distribution of charged leptons originating from W decays, to constrain 
PDFs at the percent level. In effect, W and Z production can serve as a parton-parton 
"luminosity monitor" [28]. The inferred parton-parton luminosities can then be used to 
precisely predict rates for interactions with a similar initial state as Drell-Yan production, 
such as gauge boson pair production processes. Uncertainties in the overall proton-proton 
luminosity, which is hard to measure precisely at the LHC, will cancel out in this approach. 

It is apparent from the above examples that the Tevatron and LHC physics programs 
require NNLO calculations for differential distributions in kinematic variables; knowledge of 
inclusive rates is insufficient. Although the inclusive rates for several processes, including 
Drell-Yan production of lepton pairs, are known at NNLO in QCD, until very recently 
no complete calculation of a differential quantity existed at NNLO. Such a calculation is 
quite challenging technically, and traditional methods for the computation of phase-space 
integrals cannot handle problems of this complexity. In Ref. [29] we described a powerful 
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new method of performing such calculations, and applied it to DrcU-Yan production in fixed 
target experiments. We present here in detail the computation of the rapidity distributions 
for Drell-Yan production of lepton pairs through W and {Z,j*) exchange at both the 
Tevatron and the LHC through NNLO in QCD. Although these distributions are still not 
the fully differential results needed for a Monte Carlo event generator, they allow a large 
number of the physics issues discussed above to be addressed. 

Our method extends the optical theorem to allow the tools developed for multi-loop 
computations to be applied to the calculation of differential distributions. We represent 
the rapidity constraint by an effective "propagator." This propagator is constructed so 
that when the imaginary part of the forward scattering amplitude is computed using the 
optical theorem, the "mass-shell" constraint for the "particle" described by this propagator 
is equivalent to the rapidity constraint in the phase space integration. We then use the 
methods described in Ref. [30] for computing total cross sections, keeping the fake particle 
propagator in the loop integrals, and deriving the rapidity distribution as the imaginary 
part of the forward scattering amplitude. We remark that the rapidity distribution for 
inclusive production of Higgs bosons at hadron colliders, which in the heavy top quark 
approximation is known at NLO [31], can be computed at NNLO by precisely the same 
technique; indeed, all the basic integrals encountered in the two problems are identical. 

We find that the NNLO corrections to the W and Z rapidity distributions are small 
for most values of rapidity. This is consistent with the results found in Ref. [18] for the 
inclusive cross section. However, the magnitude of the corrections can reach a few percent 
for certain invariant masses and collision energies, indicating that they are required for 
the precision desired in experimental analyses. The residual scale dependences of the 
rapidity distributions are below the percent level for all but the largest physically allowed 
rapidities. The theoretical uncertainty is therefore dominated by our imprecise knowledge 
of the PDFs. We study the effect of varying the PDF parameterization; we use several 
fits provided by both MRST and Alekhin. The different MRST sets yield results for 
the rapidity distributions that vary by Ri 1% at the LHC; the Alekhin set gives results 
that differ from those of MRST by 2 8.5% as the rapidity is varied. The anticipated 
experimental uncertainties at the LHC are sufficiently small to distinguish between such 
PDF sets. EW gauge boson production can therefore provide important information about 
the PDFs at the values of and x relevant for collider experiments. Finally, we study the 
efficacy of various approximations to the complete NNLO result. We find that the common 
approximation of including only soft gluon corrections does not accurately reproduce the 
full result for phenomenologically interesting parameter choices. 

Our paper is organized as follows. In Section H we introduce our notation. We discuss 
our method of calculation in detail in Section HL We describe the coUinear renormaliza- 
tion of the partonic cross section in Section IV. In Section V we present some analytic 
results for the partonic rapidity distributions. Numerical results for the W and Z rapidity 
distributions at both the Tevatron and the LHC are given in Section VI. We present our 
conclusions in Section VII. 
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2. Notation 



We consider the production of electroweak vector bosons V at hadron colliders, 

hi + h2^V + X, (2.1) 

where X stands for any number of additional hadrons, or partons in the perturbative 
calculation. The qiQjV coupling at the tree level is described by the interaction vertex 

Vi'j = i9vCijj''{vY + aYj5), (2.2) 
where the indices i,j denote the quark flavors: 

i,j = {u,u,d,d,...}. (2.3) 

The matrix Cij is the unity matrix when V = j,Z and is the CKM matrix when V = W. 
Numerical values for the required CKM matrix elements arc given in Section 6. 
The vector and axial coefficients for up and down type quarks are: 



2 1 

= 3> al = 0, = --, a] = 0, 

4 



= 1- - sin^ 9w, af = -l, = -I + -sin^ Ow, of = 1, 

„w _„,w _ J_ w _ w L (OA) 

The rapidity of the vector boson V is defined as 

Y^\loJ^\. (2.5) 



2 °V^-P^ 

where E and are respectively the energy and longitudinal momentum of V in the center- 
of-mass frame of the colliding hadrons. The cross section for the production of the vector 
boson can be written as the convolution of partonic hard scattering cross sections with 
hadronic parton distribution functions: 



dY ' 

da^ s--. , , „th.\, s Jho). s daY 



,^ E/ / dx,dx2f^'^Hx,)fi''^'\x2)^{x,,X2). (2.6) 

ab •'■^^ J\fTe~^ 



Here, 



■ s ■ 

my is the invariant mass of F, and S = {Pi + ^2)^ is the square of the center-of-mass 
energy of the colliding hadrons hi and h2, which carry momenta Pi and P2 respectively. 

As we will see in the next Section, it is beneficial to represent the rapidity constraint 
in a covariant form. To do so, we introduce the variable u, where 

u = (2.8) 

X2 



-4- 



In the center-of-mass frame of the two coUiding hadrons, it takes the simple Lorentz in- 
variant form 

u = P^, (2.9) 

PV ■P2 

where pi = xiPi and p2 = X2P2 arc the momenta of the incoming partons, and pv is the 
momentum of V. The partonic center-of-mass energy is s = +^'2)^ = Sx\X2- The 
partonic ^-distributions are obtained by integrating the partonic matrix elements over the 
phase space of the final-state particles with a fixed value of u: 



Here, ||A^a6^v^+x||^ denotes the square of the scattering amplitude, averaged over spins 
and colors of the colliding partons. 
The allowed values of u are 



with 



1 

z<u<-, (2.11) 



(2.12) 



S X1X2 

and 

T<z<l. (2.13) 

Inverting Eqs. (2.8) and (2.12), the arguments {xi,X2) of da^j^/dY in Eq. (2.6) are given 
by 

Xl = ^^-^= , X2 = I— . (2.14) 

^yz/u \JUZ 

The boundary values of {z^v) are only achieved for special kinematics (see Fig. 1). 
For 2; = 1, my = s, and there can be no additional partons radiated in the V boson 
production process; the kinematics is that of the Born-level process qq V. We refer to 
the limit 2; — ^ 1 as the soft limit, since any additional partons must carry little energy. The 
boundary u = z corresponds to production of a y boson along with one or more partons 
radiated collinear with incoming parton 2, with momentum (1 — z)p2- That is, inserting 
Pv = Pi + zp2 into Eq. (2.9) leads to u = z. Similarly, the boundary u = 1/z is achieved 
when the additional partonic radiation is collinear with incoming parton 1. We refer to the 
limits u ^ z and u ^ 1/z as collinear limits. 



3. Method 

We evaluate the partonic rapidity distributions of Eq. (2.10) through NNLO in perturba- 
tivc QCD. The NLO corrections have been previously evaluated in Ref. [25]. However, the 
calculation of the NNLO corrections is intractable using current techniques. We describe 
here a new method powerful enough to handle this problem. We express the rapidity con- 
straint as the mass-shell condition of a fake "particle." This permits the use of the optical 
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Figure 1: Variables {z, u) used to describe the kinematics of the vector boson rapidity distribution 
at parton level. The physical region is hatched. The point u = z — I corresponds to no additional 
radiation, or Born-level kinematics. The left edge, u — z, corresponds to radiation of partons 
collinear with incoming parton 2. The right edge, u — 1/z, corresponds to radiation coUinear 
with parton 1. The arrows show flows relevant for the convolution integrals encountered in mass 
factorization (see Section 4). 

theorem to transform the matrix elements into cut forward scattering amplitudes. We can 
then apply methods developed for multi-loop integration to evaluate these amplitudes. We 
describe in detail below the required modification of the rapidity constraint, the simplifica- 
tion of the forward scattering amplitude using integration-by-parts reduction algorithms, 
and the evaluation of the resulting master integrals. 

We begin by describing the three distinct contributions that enter at NNLO, to illus- 
trate the difficulties that arise: 

• Virtual- Virtual, which contains interferences of diagrams with only the electroweak 
boson V in the final state; 



• Real- Virtual, which contains interferences with the V boson and one additional quark 
or gluon in the final state; 
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Real-Real, which contains interferences of tree-type diagrams with V and two addi- 
tional partons in the final state. 



AAA/ ^AA^ 



The Virtual- Virtual contribution to the rapidity distribution is identical to its counter- 
part in the total cross section, which has been computed previously [18]. The new features 
of the rapidity distribution are the Real- Virtual and Real-Real components, which now 
have non-trivial kinematic constraints. Until very recently no systematic technique for their 
evaluation existed; this was the major reason for the lack of progress. However, since the 
calculation of the inclusive Drell-Yan cross section, our ability to calculate diagrams of the 
Virtual- Virtual type has progressed greatly. New algorithms for the evaluation of two-loop 
diagrams of the same [32] and more complicated topologies [33, 34] have been developed. 
It is now well understood how to organize the evaluation of generic multiloop amplitudes 
using integration by parts and Lorentz invariance reduction algorithms [34, 35, 36, 37], and 
how to compute the resulting master integrals using either the Mellin-Barnes [38, 39] or 
the differential equation method [40, 37]. Our method renders the Real- Virtual and the 
Real- Real contributions amenable to the same techniques. 

3.1 Construction of the modified forward scattering amplitude 

We follow the approach introduced in Ref. [30, 41, 31, 29]; we replace all non-trivial phase- 
space integrations by loop integrations. To accomplish this we represent all delta functions 
constraining the final-state phase space by 

^(^) = - ■ (3-1) 



2TTi \x — iO X + iO 

In the evaluation of total cross sections we only have delta functions which put the final- 
state particles on their mass shell, 

Jdllf^llJ d% 5+{p) - m}), (3.2) 

where we work in d = 4 — 2e dimensions, and 6^ includes the positive energy condition 
Ef > 0. Using Eq. (3.1), each such delta function becomes a difference of two propagators 
with opposite prescription for their imaginary part: 

Sipj-mj) - p2_l2_,Q -i^-^-)- (3-3) 

When calculating differential quantities, there are additional constraints on the phase space. 
The distribution constraints can usually be expressed through delta functions with Lorentz- 
invariant arguments that are polynomial in the momenta of the final-state particles; we 
then transform them into propagators using Eq. (3.1). 
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For the rapidity distribution of the massive boson we substitute 
'pv -pi \ PV ■P2 



u 



c.c. 



. . v^-^v (3.4) 

^Pv ■P2 J PV ■ [.Pi - UP2) - ^0 

The above substitution introduces a propagator with a scalar product in the numerator and 

a denominator hnear in the momentum of V. However, the multi-loop methods we employ 

are not sensitive to such irregularities in the form of the propagators; they only require 

that the propagator of Eq. (3.4) be polynomial in the momenta. Substituting Eqs.(3.3) 

and (3.4) into Eq. (2.10), we obtain a forward scattering amplitude with "cut" propagators 

originating from both the on-shell conditions on the final-state particles and the rapidity 

constraint. Pictorially, the three different contributions can be represented by diagrams 

similar to the following ones: 

• Virtual- Virtual; 





Real- Virtual; 



• Real-Real. 



000/ .000 
000. .000 - 



We have associated an additional "rapidity" propagator with the V boson, which we rep- 
resent by a straight line just to the right of the cut from the usual wavy (cut) V boson 
propagator. In the above diagrams, cut propagators represent differences of two complex 
conjugate terms, propagators on different sides of the cut have different prescriptions for 
their imaginary part, and the initial and final states are identical. 

The three contributions are now expressed as two-loop amplitudes in which the cuts 
denote differences of propagators with opposite ie prescriptions. These cut conditions are 
accounted for at the very end of the calculation, after using generic multi-loop methods 
to simplify the two-loop expressions. We generate the diagrams for the forward scattering 
amplitude using QGRAF [42]. We then apply the Feynman rules, introduce the rapidity 
"propagator" of Eq. (3.4), and perform color and Dirac algebra (we use conventional dimen- 
sional regularization) using FORM [43]. This generates a large number of integrals with 
cut propagators which we must evaluate. The evaluation of these integrals is discussed in 
the next Section. Our treatment of 75 in dimensional regularization follows the discussion 
in Ref. [18], to which we refer the reader for further details. 
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3.2 Reduction to Master Integrals 

An essential part of the calculation is the reduction of the integrals to a small set of 
independent master integrals using linear algebraic relations among them. This procedure 
is routinely applied in computations of virtual amplitudes, such as in the Virtual- Virtual 
contribution. Our representation of the Real- Virtual and Real-Real terms makes it possible 
to evaluate them in a similar manner. 

We first apply partial fractioning identities to reduce the denominators of the inte- 
grands to a linearly independent set. Partial fractioning is always applicable to box dia- 
grams with the introduction of the rapidity constraint. Seven independent scalar products 
can be formed from the two external momenta pi,P2 and the two loop momenta k, I (omit- 
ting the constant scalar products, pi ■ pj)- On the other hand, box topologies with the 
rapidity propagator have eight terms in the denominator. Thus the eight terms obey one 
linear relation, which can be used to perform a partial fraction decomposition whenever all 
terms entering the relation appear as denominators. This step allows us to eliminate one 
of the uncut propagators in favor of the rapidity propagator. (Whenever a cut propagator 
or the rapidity propagator does not appear in the denominator of an integral, that integral 
may be set to zero, by anticipating the delta function constraints.) Prom this procedure we 
derive 11 major topologies for the Real-Real contributions, 2 major topologies for the Real- 
Virtual contributions, and 2 for the Virtual- Virtual contributions. All non-box diagrams 
are sub-topologies of the above set. 

We obtain additional recurrence relations using integration-by-parts (IBP) identi- 
ties [35, 36]. li k,l are the loop-momenta of the two-loop integrals in the forward scattering 
amplitude, and pi , p2 are the incoming momenta, we can write 8 IBP identities of the fol- 
lowing form for each integral: 



with ri^ = k^, and = k^, l^,p'l,p^. Each integral contains some cut propagators. How- 
ever, since differentiation with respect to the loop momenta is insensitive to the prescription 
for the imaginary part of propagators, the application of IBP reduction algorithms and the 
taking of cuts commute [30] . This fact allows a straightforward derivation of reduction re- 
lations for phase-space integrals. Similarly, we can derive Lorentz invariance identities [37] 
for phase-space integrals; however, for the Drell-Yan rapidity distribution they are not 
linearly independent of the IBP relations and provide no additional information. 

To solve the system of equations formed by the IBP relations we use an algorithm 
introduced by Laporta [34]. We construct a large system of explicit IBP identities which 
we then solve using Gauss elimination. This system should include a sufficient number 
of equations to reduce all the integrals of the forward scattering amplitudes to master 
integrals. A detailed description of the algorithm can be found in the original paper of 
Laporta [34]; we have implemented a customized version in MAPLE [44] and FORM [43]. 
An important simplification of the reduction procedure in the present case is that we can 
discard all integrals in which the cut propagators, or the rapidity propagator, arc eliminated 
or appear with negative powers [30]. After performing the reductions, we obtain 5 master 




(3.5) 
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integrals for the Virtual- Virtual contributions, 5 for the Real- Virtual, and 19 for the Real- 
Real. 



3.3 Master integrals 

All Virtual- Virtual master integrals were known prior to this calculation. The only non- 
trivial one is the crossed-triangle master integral, which was calculated in Ref. [45, 46]. A 
list of all Virtual- Virtual master integrals can be found in the appendix of [30]. 
For the Real- Virtual contributions we find the following master integrals: 




where solid lines correspond to massless scalar propagators 

1 

bold solid lines correspond to massive scalar propagators 



1 

— rriy 

and dashed lines denote the rapidity propagator 



1 

PV ■ [Pl - UP2] ' 

The Real- Virtual master integrals can be evaluated by using Eq. (3.1) to reinstate the 
delta- function constraints. We then must perform a one-loop integral and a 2-particle 
phase-space integral; both are straightforward. The most complicated loop integral is a 
massless one-loop box diagram with one external leg off-shell, which is known to all orders 
in e [47]. The Real- Virtual phase-space integration is simple because the polar angle for 
the 2 — > 2 process is fixed by the rapidity constraint, leaving only a (1 — 2e)-dimensional 
azimuthal angular integration. It is thus straightforward to derive analytic expressions for 
the Real- Virtual master integrals which are valid to all orders in e. 

The Real-Real phase-space master integrals were unknown prior to this calculation. A 
few can be evaluated directly; for example, the two simplest master integrals. 



m 



(fqv(rqi(rq25 {pi + P2 - qv - qi - 92) 
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5+{ql - ml)5+{ql)5+{ql)5{qv ■ [pi - up^]) 



(3.6) 



and 




<fqv<fqid'^q25'^{pi + P2 - Qv - Ql - 92) X 



{qi + q2?5+{ql - ^1)5+ {ql)5+ {ql)5{qv ■ [pi - np2]) 
have the following hypergeometric integral representation: 



with 



r(d/2) 



and 



/C.(<5)= / dxxnx(l-x)(l-X'5)]-% 

JO 

(^ - ^)(1 - ^/ui) 



Expanding in e, we obtain 



/Co(<5) = 1 + e 



(5- l)ln(l -5) 



{5 - 1) In^ (1 - 5) 



and 



-3 



/Ci(<5) = - + - 
^ ^ 2 2 



((5-1) In (1-5) ((5 -2) Lis (5) 



+ 



6 



+ 9 



vr 



+ 0(e=^ 



(J2_i)ln(l- J) 5<5 + 2 



52 



2,5 



e2 



((52 - 1) ln2 (1 - 5) 



(3.7) 

(3.8) 

(3.9) 
(3.10) 

(3.11) 



(3.12) 



((52-1) ln(l-(5) ((52 - 2) Li2 ((5) 3,5 + 2 _ 7r2 
252 + P ^^^^"T 



52 

+ 0(6^). (3.13) 



In the expression for the scattering amplitude, some of the master integrals are mul- 
tiplied by coefficients which become singular at the phase-space boundaries. For example, 
when the /[z^] master integrals get multiplied by l/{u — z)'^^'^ or 1/(1 — uz)'^'^'^ , the matrix 
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elements become singular at u = z or u = l/z respectively. These singularities are regu- 
lated by the non-integer powers of the {u — z) and (1 — uz) prefactors in Eq. (3.8). Upon 
integrating over u and z, they generate 1/e poles which cancel against the Real- Virtual 
and Real-Real 1/e singularities. For example, 



du{u-z)-^-^' = --(^^] . (3.14) 



2e 



Since we are interested in the rapidity distribution, we do not integrate over u. We must 
therefore extract these singularities from the Real-Real master integrals. To do so, we 
factor out the leading behavior of the integral Xi in the limits u ^ z and u 1/z, keeping 
the exact e-dependence: 

Xi{z, u) = {u- z)"^-"^(l - uz)'^-f^'J^i{z, u). (3.15) 

The integers m,n are characteristic to each master integral, while a = /? = 2 for all 
Real-Real phase-space integrals. The functions Ti are smooth and non-zero at u = z and 

u = 1/z, and can be calculated as a series in e. In the non-singular regions of phase 
space we need only calculate the first few terms in the e expansion, up to the order where 
polylogarithms of rank 2 appear. However, at u = z and u = 1/z additional 1/e coefficients 
may be generated, and at u = z = 1 additional 1/e^ poles may appear. These require an 
e expansion of J^i up to a transcendentality of rank 3 or 4. We therefore split the master 
integrals into four different terms: 

Mi = ;ff + ;f^=°"(") + -Y™"^^^ + Xl^'^''^. (3.16) 

Here, 

= {u- z)"^-"^(l - -uz)"-^^J^j(l, 1) (3.17) 
is potentially singular at the limits u = z, u = 1/z and u = z = 1; 

can only become singular at u = z; 

^coui-) = [--zj (1 - uzT-^' [Ti{z, 1/z) - Tiil, 1)] (3.19) 
can only become singular at u = 1/ z\ and 

is smooth in all singular limits. We extract the explicit 1/e terms from the coll and soft 
terms by replacing the u variable with 
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where <y <1. We then apply identities of the form 



X 



-l+e 



=-six)+Y: 



(3.22) 



+ 



for a; = y, 1 — y and \ — z. The advantage of using the variable y instead of u is that y 
separates the singularities at u = z and u = l/z, which overlap when z = \. [At next-to- 
Icading-ordcr, and for the Real- Virtual 2-particle phase space at NNLO, the variable y is 
related to the 2^2 partonic center-of-mass scattering angle 6** by y = (1 + cos^*)/2.] 

Although a deeper expansion in e is required for the master integrals in the collinear 
and soft regions, the calculation is simplified since in the collinear regions the result has a 
non-trivial dependence on only the variable z\ in the soft region, the Ti have no dependence 
on either u or z. For example, while it is difficult to expand K^yiS) to higher orders in e for 
generic (5, in the soft region, 5 — > 0, it can be computed in terms of Gamma functions in 
closed form: 

r(l + - e)r(l - e) 



/C.(0) 



r(2 + 1/ - 2e) 



(3.23) 



3.4 The differential equation method 



The two Real-Real master integrals of the previous subsection were calculated by deriving a 
simple hypergeometric integral representation starting from their definition as phase-space 

integrals. However, this is not practical for most master integrals. In more complicated 
cases we resort to the method of differential equations. This method was developed for 
loop integrals [40, 37]; however the representation of Eq. (3.1) for delta function constraints 
allows its application to phase-space integrations in a straightforward manner [30]. We 
consider the following master integral as an example: 



J{z, u) = J d'^qvd'^q\d'^q2 S'^{pi +P2-qv-qi- ^2) x 
S+iql - ml)S-^{ql)5+iqi)d{qv ■ [pi - UP2])- 



1 



(qi + qv -piy 

After applying the transformation of Eq. (3.1), this integral becomes 



(3.24) 



J{z,u) = J 



d'^k dH 



1 



fc2 



m 



V 



1 



-1 



.{l+Pi +P2Y 
where k = —qv, I = —qv — qi, SLud we denote 



k- {pi- up2)]c (l+Pi) 



1 



^2' 



'1" 






_x_ 


c 


2TTi \ 



2m \x — ■iO X + iO^ 
We can now differentiate J{z, u) with respect to z and u, obtaining 



(3.25) 



(3.26) 



= J d'k dH 



1 



(A;2 — myY 



1 



{k-lY_ 
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du 



{I+P1+P2Y 



1 



-1 

k-{pi- UP2) 
1 



(3.27) 



k"^ — my 



{k-l)\ 
-1 



k-p2 

{I+PI+ P2f\ c[{k-{pi- Up2)f\ ,{l+Pl?' 



(3.28) 



We have set s = (pi +P2) = 1 in these expressions. Neither integral on the right-hand 
side of Eqs. (3.27) and (3.28) is a master integral. However, using IBP we can reduce them 
to the master integrals J, /[O], and J[l], using the reduction algorithm of Section 3.2. We 
then obtain a system of two partial differential equations which determines the functional 
dependence of J on the two kinematic variables z, u: 



dj{z,u) 
dz 



dj{z,u) 
du 



2e ^, , il-2e)u\il + ^zu + Az)e-l-zu-2z\^,,, 



u — z 

(1 - 2e)(2 - 3e)u 



2ez{u — z){l — uz) 



2ez{u — z){l — uz) 
I[l]{z,u), 



(3.29) 



2e , (1 - 2e) [(7 + Az - 3zu)e - 3 - 2z + zu] , , 
Ji^^ ^) + ^ n.r.. ^^[0](^, u) 



+ 



u — z 

(1 - 2e)(2 - 3e) 



2e(u — z){l — uz) 



2e{u — z){l — uz) 
The general solution of Eq. (3.29) is 

1 2 



I[l]iz,u). 



(3.30) 



J{z,u) 



2d-2 



{u - z) 



-2e 



I J' dzi {u - z^f P{zi,u) + fiu) + c| , (3. 



31) 



where (Od_i/2'^^^)^ f3{z, u) is the inhomogeneous part of the differential equation in Eq. (3.29). 
We can evaluate the integral in Eq. (3.31) as a series in e after we rewrite P using the ex- 
pressions for from subsection 3.3. We obtain 

n = J' dzi {u - zif P{zi,u) = + Ao{z, u) + 0(e), (3.32) 



with 



and 



Ao{z,u) 



Ai {z, u) = - Mr) - ln(r + t)] (3.33) 



[ln(r) + ln(2)] ln(l + r^) + ^ ln(r) [21n(r) + ln(2) - 4] 
~ \n{t) ln(r) + ^ [ln(2) + 4 - ln(r) + ln{r^ + 1)] ln(r + *) - ^ ln^(r + t) 



+ ^ln(t)ln(r + t) -Li2 



+Li2 



2r' 



:Lio 



(r + t)r 
r2 + 1 
r-f 
2r 



+ -U2 
2 



:Lio 



t 

r 

{r — t)r 
r2 + 1 



+ -Li2 
2 



r-t 



(3.34) 
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We have introduced the notation r = ^/u and t = \[z. Substituting the solution of 
Eq. (3.31) into the differential equation of Eq. (3.30), we derive a differential equation for 
j{u) which we can again solve order by order in e. We find 

/(^) = Mf) +/o(^) + 0(e), (3.35) 



with 



and 



/i(^x) = ^ln(l + r2)-In(r) (3.36) 



hiu) = Li2 [-r^] + Li2 ^^^^ - Li2 [r^] + 1 In^r"^ + 1) - ln2(r) + 41n(r) 

[ln{r^ + 1) + 2 ln(r)] In(2) + [ln(r) - 2] ln{r^ + 1). (3.37) 

Finally, we must determine the constant of integration C. In principle, this requires an 
explicit calculation at a specific kinematic point J^{zo, uq)- However, in many cases we can 
extract the constant of integration by comparing to the asymptotic behavior of all rapidity 
phase-space integrals &t u = z = 1, which is identical to that of the basic master integral 
/[O]: 

lim PS(z, u) = c{u - zy'-^'{l - uzy-'^^ (3.38) 

2,M— »1 

The e power of the u — z and 1 — uz factors is determined by the number of dimensions, 
cZ = 4 — 2e; adding more propagators to the basic master integral I[0] can only alter the 
integers n,m of the asymptotic scaling. We note that the presence of the constant of 
integration C in Eq. (3.31) violates the scaling of Eq. (3.38). We can therefore evaluate C 
by requiring that all the terms in Eq. (3.31) that violate Eq. (3.38) in the limit z — 1, u — 1 
cancel. We obtain 



2- 

■ + 0(e). (3.39) 



2 

There are master integrals for which the solution of the homogeneous differential equation 
gives a scaling at n = z = 1 which is consistent with Eq. (3.38) for arbitrary values of the 
constant C. For these master integrals, we must determine C by performing an explicit 
evaluation in the vicinity of this kinematic point. 

As discussed previously, we often need to calculate master integrals in their soft or 
collinear limits to higher orders in e. For example, the integral is typically divided by an 
explicit (n — z) factor in the matrix elements, requiring an e-expansion in its collinear limit 
u ^ z which includes the order e term. We could extend the outlined calculation of J' for 
generic z,u to include the 0{e) term and then take the limit u ^ z. However, this would 
involve expressing the result for generic z, u through generalized poly logarithms of rank 3 
with two variables; taking the u ^ z limit would collapse them to rank 3 polylogarithms 
with only the argument z. We can avoid the two-variable rank 3 polylogarithms by solving 
the differential equations directly in the u ^ z limit. We express the z-dependent term in 
the general solution of Eq. (3.31) in the form 



7^ 



r dzi{u-zif'P{zi,u), (3.40) 

J z 
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and perform the change of variables 

zi = z+{u- z)\. (3.41) 

Next we expand the integrand in u — z and keep only the leading term. Only the coll{z) 
limits of the boundary integrals /[i/] are required, and as explained above those are known 
to all orders in e. We can then expand Eq. (3.40) in e; the resulting integration over A 
involves polylogarithms with a single argument z, and can be performed straightforwardly. 
The computation of lim„^^/(tt) proceeds as before, utilizing equivalent expansions in 
u — z. Finally, the constant C is determined by matching to the asymptotic behavior, 
(u-z)-2^(l-z2)-2^ 

An important check of our results for the master integrals is provided by integrating 
them over the rapidity variable u. The master integrals also enter the NNLO corrections 
to the rapidity distribution for Higgs boson production at hadron colliders via gluon-gluon 
fusion, computed in the heavy top quark approximation. Hence the integrated master 
integrals can be expressed in terms of the master integrals appearing in the evaluation of 
the Higgs boson total cross section. We have verified that all rapidity-distribution master 
integrals are consistent with the results of Ref. [30]. The analytic expressions for the 
master integrals are too lengthy to present here. They can be obtained from the authors 
by request. 

4. Renormalization and mass factorization 

The partonic cross sections of Eq. (2.6), after combining the real and virtual contributions 
up to 0{a1), contain and 1/e poles arising from both ultraviolet and initial-state 
collinear singularities. We remove the UV singularities through renormalization in the MS 
scheme, and absorb the initial-state singularities into the PDFs using the MS factorization 
scheme. First we expand the cross section in the strong coupling constant, 



The bare strong coupling a'g is related to the running strong coupling constant ag = ag (/i) 
in the MS scheme via 



a'g {Any e-^^ = a, /x^ 

with 



l_^^ + 0(a2) 
TT e 



(4.2) 



Here is the number of light quark flavors, and fi = jiR = fip is the combined renormal- 
ization and factorization scale. At the end of the calculation, we restore the dependence 
on alone, with the aid of the renormalization group equation. Substituting Eq. (4.2) 
into Eq. (4.1) and collecting with respect to ag, gives the coefficients of the renormalized 
expansion, 



dY dY KttJ dY \nJ dY 



+ ^ -7^ + ^ -J^ + 0{{agn (4.4) 
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in terms of the bare ones in Eq. (4.1). 

Similarly, to remove the initial-state singularities, we rewrite the hadronic cross section 
of Eq. (2.6) using infrared-finite partonic cross sections: 



ab 



'^"'^ Et t dx,dx,h^^\x,)fi^^\x/^{x,,x,). (4.5) 



The renormalized parton distribution functions fa^^ are related to the "bare" ones /^^^ by 

/f^ = /f^®r„,. (4.6) 

We have introduced the convolution integral 

(/ ® g) (x) = [ dydz f{y)g{z)5{x - yz), (4.7) 
Jo 

and we implicitly sum over repeated parton indices. The functions Faj, are given in the MS 
scheme by 

Tabix) = SabS{l - X) - 

TT e 

+ (?)' {2? [(^« ® ^i?') + H - ^^""<^'} + °<°')- 

where the Altarelli-Parisi kernels P^^^ can be found in Refs. [48]. Substituting Eq. (4.6) 
into Eq. (4.5) and comparing with Eq. (2.6) we find 

^(^,u)=r dy, r dy,TUyi)^-^(^,^-^)rM. (4.9) 
dY dY \y1y2 y2 J 

The convolution integrals follow contours in the {z,u) plane, as shown in Fig. 1. The yi 
integration, holding y2 fixed, sweeps out a flow such as the one marked "left collinear"; 
whereas the ?/2 integration sweeps along a "right collinear" line. The lower limits of the 
integration over the t/j correspond to the point {z,u) = (^^, ^u) striking one of the two 
boundaries, u = z or u = 1/z. We solve Eq. (4.9) for the finite partonic cross sections 
duab/dY recursively, order by order in the expansion. 

At this point it is straightforward to derive the finite partonic cross sections. We 
outline below the salient features of the calculation. All cross sections referred to in the 
formulas below are finite; we henceforth drop the superscript tilde when referring to them. 
We will also drop ^^d/dY" to make the formulas more compact. 

• To O (og) , at least one of the two Tab factors, or daab/dY, on the right-hand side of 
Eq. (4.9) has a delta function containing the convolution variable. If neither Tab factor 
contains a delta function, then only the LO cross section enters, with u = z = 1. This 
forces both yi and y2 to be set to their lower endpoints, sjzju and ^Jzu respectively. 
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so no integral needs to be done. Apart from this case, the double integral in Eq. (4.9) 
reduces to a single integral of one of the following two forms: a "right" convolution, 



or a "left" convolution, 



{z,u)= [' dxaab(-,-)P^\x), (4.10) 

\/uz ^ 



{z,u) 



J ^dxaab{^,xu)pl;!\x). (4.11) 



Using the behavior of the partonic cross sections under inversion of rapidity, aab {z, ^) 
aba {z, u) , it is simple to show that 



'^ab<^pi:^ 



1 

Z, - 
U 



P 



in) 



be 



(^ba 



{z,u) . 



(4.12) 



We need only consider right convolutions; we can obtain left convolutions by inverting 
the variable u. 



• The convolutions required to obtain a finite NLO cross section are of the form a 



P^ ■ The LO cross section has the following form: 

afj ^6{l-z){5{y)+6{l-y)}, 



ab 



(4.13) 



where we have used the variable y defined in Eq. (3.21). Substituting a^*^^ into the 
convolution formula in Eq. (4.10), we find that the resulting 5 (l — |) removes the 
integration, leaving only the product of {z) with the remainder of a^^h We 
note that this remainder contains either 6 (y) or 6 {1 — y). To put it another way, in 
Eq. (4.10), since a^^ {z/x,u/x) requires z/x = u/x = 1, the terms generated all have 
u = z, corresponding to y = 0. The S {1 — y) term only contributes in the limit of 
Born kinematics, u = z = 1. 

• There arc three distinct types of convolutions needed in the NNLO cross section: 



^ib ' 



,(1) 
be ' 



(0) ^ p(0) 



ab 



be 



. The first of these is simple; as in the 



NLO cross section, the ^ (l ~ f ) from the Born cross section removes the convolution 
integral, and all the terms generated have u = z. We discuss the remaining cases 
below in some detail. 

— We solve the second type iteratively. The o'^J ^Pj^^^ piece was already computed 
to obtain the NLO cross section. It contains either S (y) or 5 {1 — y), as noted 
above. It may also contain distributions in 1 — z. It is simple to show that 
when performing the second convolution integral using Eq. (4.10), 5 (1 — y) — ^ 
S {x — ^/uz) (and again u = z), removing the integration. In the S (y) terms, it 
is convenient to treat plus distributions as follows: for distributions of 1 — z, we 
set 



'In" (1 - z) 
1-z 



(1 



-l+€ 



(4.14) 
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where the vertical bar indicates that we should take the appropriate term in the 
e expansion defined in Eq. (3.22). For distributions of 1 — x arising from the 
splitting function, we use 

(4.15) 

a" 

where we now must take the O (a°) term. The most complicated integral we 
must evaluate, which contains plus distributions in both 1 — z and 1 — x, becomes 

h = 5 (y) £ dx fiz/x) (l - ^) (1 - x)-'+^' , (4.16) 

where f{z/x) is finite in all kinematic limits and we have used the delta function 
to simplify the lower limit of integration. Performing the variable change q = 
{x — z)/{l — z), we obtain 

h = 5 (y) /' dq [q{l - z) + z]^'' f{z/x[q]) (1 - z)-'+<'+^^ q-'+' (1 - 9)-^+"^ , 
Jo 

(4.17) 

where x[q] = q{l — z) + z. We can extract the distributions in 1 — 2; by using 
the expansion in Eq. (3.22). We must also interpret the q and 1 — q factors as 
distributions; we set 

1 

q-^+^=U{q)+J2- 
n 

and utilize a similar expansion for ^- — q. We can now expand the integrand 
to O (0*^6™) . Performing the required integrations, we obtain the result for this 
convolution in terms of poly logarithms of rank 2 and 3 in the variable z. 

— To obtain convolutions of the form , we first return to the form of the 
NLO cross section before expansion in e, which is 

oc y-i-^ (1 - y)-^-' (1 - z)-^-^^ + . . . . (4.19) 

The ellipsis denotes terms of the form CF^ah^-^bc ' which are needed for an infrared 
finite NLO cross section; the convolution of these with pj^^ has already been 
discussed, and we ignore them here. We have presented the qq cross section; the 
qg NLO result differs only in the exponents of j/, 1 — y, and 1 — z which appear, 
and the required convolutions proceed similarly to those we now discuss. We 
again consider the case where the splitting function contains a plus distribution 
in 1 — X. We rewrite this term using Eq. (4.15). The integral we must evaluate 
becomes 

l2= f y;i-^(l-yp)-^-^(l--)"'"''(l-x)-^+"% (4.20) 



\-X 



(1-x) 



-1+ae 



In" 



(4.18) 
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where 

_ x{u-z) 

{x-z){x+u) ^ > 

and f{§,yp) is finite in all kinematic limits. Performing the variable change 
q= {x — s/uz) / (1 — V^)) integral becomes 

h = C dq fiz, y; x[q]) g{z, y; x[q],e) y"^"^ (1 - y)-i-(i-) (1 - 
JO 

xg-i-^(l-g)-^+"% (4.22) 

where x[q] = q{l — \/uz) + ^/uz. We have absorbed terms which are finite 
in all limits into the function g. We extract the singularities in y, 1 — y, and 
1 — z using the expansion of Eq. (3.22); we again interpret the q and 1 — q 
factors as distributions, and expand them as in Eq. (4.18). We can now expand 
the integrand in both a and e. To obtain the contribution to the NNLO cross 
section, we take the O (a*^) term, and expand it in e up to and including the 
O [e^) piece. The resulting integrals are straightforward to evaluate, and again 
give polylogarithms of ranks 2 and 3. The rank 3 polylogarithms only appear 
in the d{y) terms, and are functions of z only. 

After performing both the UV renormalization and the collinear subtractions discussed 
above, we obtain finite partonic cross sections. 

5. Pcirtonic cross sections 

The basic quantities we compute, (fa^^^^^^°"^/dM/dY, include the probability for the 
vector boson V to decay into a pair of leptons, e.g. Z — or — I'^ui, and are 
differential in both rapidity Y and di-lepton invariant mass M. Wc shall present our results 
in a format which is normalized properly for virtual photon production, 7* (see 
Eq. (6.2) below). For W and Z production, as well as for 7-Z interference in the 
channel, we introduce additional normalization factors A^^, where: 

ATT = 1; 

z 3 TzBf 



lesf^^cf^OQED Mz (M2-M|)2 + r|M|' 



vjvf M'^{M'^-Mf 
Ssl^cl, (M2 - M|)2 + TlMz 



2 ' 



4s2^,aQED Mw (M2-M2,)2 + r2^M2,- 

We have used the notations Vz and Tw for the total widths of the Z and W, Mz and Mw 
for their masses, and Bf and BJ^ for their branching fractions into leptons. The leptonic 
vector couplings appearing in ATT^ are given by 

v'J = -l, vf = -l + Asl., (5.2) 
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and sw and cw represent the sine and cosine of the weak mixing angle, respectively. 

Finally, we require the luminosity functions LJj {xi,X2) that enter the hadronic rapid- 
ity distribution. These functions contain the PDFs for the partons and appropriate 
combinations of the electroweak couplings to V. We follow closely the notation of Ref. [18]. 
We first introduce the following 2n/ x 2nf matrices: 



C'^± iqk,Qi) 
Cw± {qk,qi) 

c^± iqk,qi) 



lifqk = qi 

otherwise ' 

1 if % = « 

otherwise ' 

1^9*5; P if e,, +eg; = ±1 
otherwise 

if eg, = ±l + eg, 
otherwise 

I^MiP if eg, +eg( = ^1 







otherwise 



(5.3) 



Here is an element of either of the following j^-dimensional vectors: Q = {u,d,s,c,b}, 
Q = {u,d,s,c,b}. In Eq. (5.3), Cg^. denotes the electric charge of the element, and V^j,g, 
indicates the appropriate CKM matrix element. Using these matrices, we can write the 
luminosity functions as follows: 



L 



BC 



^AB,vec 



V 
19 

V 
9Q 

-V 

JQ2 



L 



V 
D2 



L. 



V 

CD,vec 



xi,X2) = ^ Cv {qi, qj) (v^'"^ + a^'"^^ qi (xi) qj (^2) ; 



i,j€Q,C 



xi,X2) = ^ ^ C{/ {qk, qi) (v^''^ + a^'"^^ qt (xi) qi (^2) ; 

i&Q,Q k,leQ 

xi,X2) = ^ ^ [Cy iqi,qk) + Cv {qi,qk) (v^'^ + a['^^ qi (xi) qi {x2) ; 

ieQ,Q k€Q,Q 

XI X ilk: qk) v(v\ qi (xi) qi {X2) ; 

x\,X2) = X qj) (^]^'^ + a]^'^) qi (a^i) 5^ (3^2) ; 

xx,X2) = L^g {X2,xi) ; 



Xi,X2) 



Xi,X2 



Xi,X2) = X ili^ ^k) {vj''^ + qi {xi) qj (^2) ; 

i,j&Q,Q k£Q,Q 

xi,X2) = X ^ ilvik) iv]''^ + aj'^) qi {xi) qj {X2) ; 



i,j&Q,Q k€Q,C 

Xl,X2) = Yl ^ ili: 1i) '^Y'^J 1i (^1) I3 (2^2) ; 



i,j&Q,Q k&Q,Q 
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LcEi {X1,X2) 



i,jeQ,Q keQ,Q 

Y ili' I3) {^Y^ + (h'^) ii (^1) I3 (^2) ; 



i,j£Q,Q 
V 



LcF {xi,X2) = Y Yj I3) {"T'^ + 1i (^1) 1i (^2) ; 

L^g {xi,X2) = Y ^3) (^i^'^ + Oi^'^) 9 {Xl) 9 {X2) • 



(5.4) 



In this formula, a function such as qi (xi) denotes the appropriate parton distribution 
function. The label V takes the values 7, Z, and W^. The electroweak couplings vf and 
af are given in Eq. (2.4). To obtain the 7-Z interference luminosity functions, we must 
use V = J and substitute vj''^ — vjvf, vjvj — 5 (^ivf + i^J^f^ ■ 

The final ingredients required are the partonic hard cross sections for the channels cor- 
responding to the luminosity functions (5.4), (io"]^/(iy(z, u) for ij G {NS, B^, BC, . . . , gg}. 
We have obtained analytic expressions for these functions; however, they are quite lengthy, 
so we refrain from giving them here. A MAPLE file containing the functions is available 
from the authors by request. They have also been implemented in C-|— 1-, as part of a 
numerical program computing the hadronic rapidity distribution. The bulk of the ana- 
lytical complexity stems from the "hard" region, away from the boundaries, z < 1 and 
z <u <l/z (or < y < 1). 

The hard functions contain polylogarithms of rank 2, Li2{Ai{z,u)), and there are a 
large number of possible ways the arguments Ai can depend on the underlying variables 
z, u. In most cases, the arguments are rational functions of i = ^/z and r = ^/u, as 
in the case of the sample integral J{z,u) presented in Eqs. (3.34) and (3.37). In four 
cases, though, we have to introduce functions in which the polylogarithmic arguments are 
significantly more complicated. The four functions of this type, J3, J27, J21 and J2, are 
given by: 



J'iiz^u) 



1 



-Re 



-^ln^(z/u) - ^ln(l + t()( \n{z/u) + 21n(l +tr) 



+ Li2 
— Li2 



-2 In 

2t(l + u) 

di — r 
2z{l + u] 
r{d\ — r) 



di-r- 2t{l + u) 
di — r 

V di+r 

Li2( + < 

r{di + r) 



2 In 



di + r + 2t{l + u) 
di+r 



(5.5) 



where di = y^u + 4z{l + u); 



J27iz,u) = - 



2rd^ 



-Re 



In 



r2 — ri + rir2 — 3 — 2u 
r2-ri- rir2 + 3 + 2m 
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+ ln 
+3 In 

where ri = vT-P~u, r2 = y/5 + Au; 



di - 2tri + rir - ndi +2t + 2tu + r 
di — 2tri + rir + ridi — 2t — 2tu — r 
{r - di + 2tri){l + r2 - 2ri) \ 
(r + di-2iri)(l-r2 + 2ri)/ 



(5.6) 



•^21 It) = Re< 



2;(1 + u)xi 



ln{tal) In 



- 1 
1 - trat 



+ lii(ta^ ) In 



— 1 

1 — traT 



- Li2(l - a^-) + Li2(l - ira|) + Li2(l - a^) - Li2(l - tro^ 



(l + z)(l -tr) 



(1 - 2;)2(1 + u)r(r + t) 



Inz In 



1 + z 



(5.7) 



where 



2u 



1 + u J z 



2u 



1 



1 + z' 
2u 
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where 0(x) is the Heaviside function, 
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After the use of poly logarithmic identities, the set of arguments of the remaining 
polylogarithms, Li2(A), can be reduced, if desired, to 

J 1 + z 1 — u 1 — u t t r — t 1 — trl — rr — 1 

1 — t r(l — t) u — z 1 — uz 1 — tr r{r — t) u — 1 1 — n 



r + 1' r + 1 '1 + n'l + u'l + u' l + u ' r{r + t)' l + tr \ ' ^^'"^^^ 
The arguments of the logarithms that appear, In(Bj), are drawn from a simpler set, 

Bi G {z,l- z, 1+ t, 1 + z,u,r -1, r + 1,1 + u,r -t,r + t,l-tr,l + tr}, (5.13) 

but since they can appear in pairs, there are still quite a few terms of the form In(Sj) ln{Bj). 

As mentioned above, rank 3 polylogarithms of a single variable, z, are generated in the 
coUinear regions, u = z {5{y) terms) and u = 1/z (5(1 — y) terms). These coUinear terms 
have a similar form to the NNLO total cross section, integrated over rapidity [18, 19]. We 
can write the functions appearing, Li3(aj), in terms of 

r ^ ^ ^ ^l + zl-zl-z2z 

Oj G < z, —z, I — z, —1 — z,l — z , 



2 ' 2 ' 1 + z' 1 + z 
1 + z' 



(5.14) 



2'2(l + z)' 2z'2(l + z)' z 
The rank 2 polylogarithms appearing in the coUinear terms, Li2{bi), have arguments 

bi e |z, -1 - , -| , -1 , -l±i| , (5.15) 

while the arguments of the logarithms, In(ci), are 

Cie{z,l-z,l + z,2 + z,l + 2z}. (5.16) 

The hard functions have integrable, logarithmic singularities in the soft limit z ^ 1 
and the coUinear limits u — > z and u 1/z. However, the complexity of the analytical 
formulas is such that many of the individual terms in the hard functions have much more 
severe singularities in these limits, e.g. several powers of 1/(1 — z) as z ^ 1. These spurious 
singularities lead to unacceptable roundoff error. For this reason we construct patching 
functions, which are used instead of the full functions in thin strips near the singular 
regions. The patching functions are typically constructed by taking the appropriate limits 
analytically. Fig. 2 shows the regions in the (z, u) plane which have to be patched. In 
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addition to the soft and collinear regions, there are two other types of regions where the 
singularities are completely unphysical. For z = [2u/(l + u)]^, the variable x\ in Eq. (5.8) 
vanishes, leading to a singularity in functions containing J2\{z^u). There is an equivalent 
singularity at z = [2/(1 + u)]'^ in functions containing J2i{z, 1/u). In this pair of strips, 
the true function is smooth enough that an analytic patch is not necessary; instead, when 
the point (z, u) lies in the strip, we replace its value by the average of two nearby values 
on either edge of the strip [49]. Finally, the limit n — > 1 is singular, as indicated by the 
presence of (r — 1) in the set Bi in Eq. (5.13); there are spurious power-law singularities as 
well in this limit. 
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Figure 2: Regions in the {z,u) plane for which the hard functions have to be patched, because 
of singular behavior. Besides the soft limit z ^ 1, and the left and right collinear edges, there are 
spurious singularities as u ^ 1, and as z ^ [2'u/(l + u)]^ and z [2/(1 + w)]^. 

The expansion of the hard functions in a series about z = 1 can be carried out to 
very high order, and produces an approximation to the integrand which is free of spurious 
singularities. Working to order (1 — z)^^ results in an expression whose accuracy is com- 
pletely adequate for predictions for typical fixed-target kinematics [29] and for W and Z 
production at the Tevatron. However, for the case of W and Z production at the LHC, the 
small value of r = M^/5 4 x 10"^ means that values of z ~ 0.01 are actually relevant in 
the numerical integration. We have therefore used the exact, unexpanded, representations 
of the hard functions (plus patches) in order to get sufficient accuracy for the case of the 
LHC. 

6. Numerical results 

In this section we present numerical results for the W and Z rapidity distributions at 



Born kinematics 
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both the Tevatron and the LHC. We use the foUowing parameters: Mz = 91.1876 GeV, 
Vz = 2.4952 GeV, Bf = 0.03363, Mw = 80.426 GeV, = 2.118 GeV, and B]^ = 
0.1082. We use the Z-pole value of aQED(Af^) = 1/128 for the fine structure constant, 
and set sin^Ow = 0.23143, the effective mixing angle measured in Z pole asymmetries at 
LEP and SLC [50]. We expect that these choices account for the bulk of the factorizable 
electroweak radiative corrections, which dominate for nearly resonant production of W 
and Z bosons. A more accurate description would require a consistent accounting of the 
electroweak corrections [11]. 

We also need the following values of the CKM matrix elements to compute the W 
cross section: 



The absolute values of tlie other matrix elements are obtained by requiring unitarity of the 
CKM matrix. Because the collider center-of-mass energy is large, it is possible in principle 
to produce top quarks in association with the W or Z; however, since these processes can 
be distinguished experimentally, we exclude them from consideration. We also omit top 
quarks from the virtual corrections, and set the number of light (masslcss) quark flavors 
n/ to 5 in all numerical results in this paper. At one loop, the partonic subprocesses 
qq — Zg and qg — > Zq include triangle graphs, weighted by the axial couplings for the 
quarks circulating in the loop. For massless quarks, these contributions cancel generation 
by generation. The effect of a finite top quark mass on the t — b contribution has been 
studied previously, and found to be negligibly small [18, 51], so wc omit it here. 

In the previous sections we discussed how the rapidity distributions of electroweak 
bosons in partonic collisions can be computed. To obtain results for hadronic collisions, we 
must convolute the partonic differential cross sections with parton distribution functions 
which describe the probability of finding a parton with a given momentum fraction inside 
the hadron. The corresponding formula reads: 



where daYj / dY is the partonic cross section, N'^ is the normalization factor for the boson 
V, and Lfj {xi,X2) is the corresponding luminosity function; these were discussed in the 
previous section. There are three observable cross sections: production of a W~^; production 
of a W~; and neutral current production of a lepton pair , which receives contributions 
from both 7 and Z exchange as well as from 7-Z interference. 

It is convenient to change the integration variables in the above formula and express 
the integration over xi and X2 through the partonic variables z and y. Consider the case 
of negative rapidity Y; the results for F > can be obtained by substituting Y —Y 
in the formulae below. For Y < 0, using the relations (2.8), (2.12), (2.14) and (3.21), the 
integration over xi and X2 in Eq. (6.2) can be rewritten as: 



|Kd|= 0.975, IK 



0.222, \V,d\ = 0.222, 



0.974. 



(6.1) 




(6.2) 




1 



1 
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dz j dy Fij{z,y), (6.3) 
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(6.4) 



This representation is convenient for numerical integration. 

We now present results for the W and Z rapidity distributions. For the NNLO calcu- 
lations, wc use the corresponding set of MRST parton distribution functions. The MRST 
code contains four different sets of PDFs. As mentioned in the inroduction, the complete 
NNLO evolution kernels needed for a consistent extraction of PDFs at NNLO are not yet 
known. The MRST program contains both the fastest and slowest possible perturbative 
evolutions, based upon the known moments of the required DGLAP equations. A third 
set allows an evolution between these extremes. Finally, a fourth PDF set which seems 
preferred by large Et jet production at the Tcvatron is included. Unless stated otherwise, 
we use mode 1 of the MRST NNLO PDF code, which corresponds to the intermediate rate 
of evolution. 

For the most part, we present double-differential cross sections, including the decay to 
leptons, 

^2^V'^leptons 

dMdY ■ ^'-'^ 
For the case of on-shell vector bosons, these are evaluated at the resonance peak, M = Mw 
or M = Mz- Of course any experiment will integrate over the resonance profile. If this 
integral is performed in the narrow-resonance approximation, and if the 7 exchange and 
7-Z interference terms are neglected in the case of the Z, the result is 



da^ 
dY 



V 



2^^ dMdY 



(6.6) 



M=Mv 



The narrow-resonance conversion factor, TrTvf^, numerically evaluates to 3.919 GeV for 
the Z boson, and 3.327 GeV for the W. One can further integrate Eq. (6.6) over the 
rapidity Y to obtain the theoretical prediction for the "total cross section times branching 
ratio," x Bj^ ■ Our total cross section results for the MRST PDFs, for example, agree 
with results obtained using the numerical program of Ref. [18], after we omit b quarks 
from the initial state [52, 53]. (Wc note that Eqs. (B.13) and (B.16) in the article in 
Ref. [18] are missing a factor of Tf = i, and the "103" at the end of Eq. (B.ll) should 
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have an x multiplying it. Also the normalization of the W cross section in Eqs. (A. 3) 
and (A. 11) should be a factor of 2 larger. All these factors arc properly included in the 
numerical program [18].) Our program is also capable of integrating over a range of di- 
lepton invariant masses, without making the narrow-resonance approximation, and we shall 
present one such plot below. 

We first present, in Fig. 3, the rapidity distribution for a Z boson produced on-shell 
at the LHC. The LO, NLO, and NNLO results have been included. We have equated the 
renormalization and factorization scales, and have varied them in the range Mz/2 < n < 
2Mz- At LO the scale variation is large, ranging from 30% at central rapidities to 25% at 
y Ri 3. This is reduced to 6% at NLO for all rapidities. At NNLO, the prediction for 
central rapidities stabilizes dramatically; the scale variation is ~ 0.6%. This increases to 
1% at y ~ 3, and 3% at y ~ 4. However, it seems that for y < 3, the rapidity values 
accessible in LHC experiments, the residual scale dependence is no longer a significant 
theoretical uncertainty when the NNLO corrections are included. 

The magnitude of the higher-order corrections exhibits a pattern similar to that of 
the scale variation. The NLO corrections significantly increase the LO prediction; the LO 
result is increased by 30% at central rapidities, and by 15% for larger rapidity values. 
They also change the shape of the distribution, creating a broad peak at central rapidities, 
as is visible in Fig. 3. The results stabilize completely at NNLO. The NNLO corrections 
decrease the NLO result by only 1-2%, and do not affect the shape of the distribution. 

For most of the plots in the paper, in order to estimate the uncertainties in the NNLO 
predictions we shall continue to set /if = fJ'R = fJ' and vary the common scale /j, from M/2 
to 2M. However, it is useful to consider a broader range of scale variations, for at least 
one kinematic configuration. In Fig. 4 we study dependence on //j? and /ir in more detail 
for the case of on-shell Z boson production at the LHC, at the precisely central rapidity 
point y = 0. For each order in perturbation theory (LO, NLO, NNLO), using the MRST 
PDF sets we plot three curves, corresponding to 

• Common variation of the renormalization and factorization scales, fJ,F = fJ,R = fj,, but 
over a larger range of fi, M/5 < /x < 5M (solid curves); 

• variation of the factorization scale alone, setting /xr = Mz (dashed curves); 

• variation of the renormalization scale alone, setting /ip = Mz (dotted curves). 

Because the LO result is independent of agifiR), the third curve is trivially constant at LO, 
and the former two LO curves lie on top of each other. We can see from Fig. 4 that the 
tiny NNLO scale variation in Fig. 3 is not peculiar to the range M/2 < /x < 2M used there. 
Even extending the range to M/5 < fi < 5M, for a common variation the bandwidth only 
enlarges from 0.5% to 1.2%. Over this same range, holding fip fixed and varying /ir also 
produces a quite small range of values, less than 0.5%. The largest variations arc found 
by holding hr fixed and varying /ip. These variations are still only of order 0.7% over 
the range M/2 < fi < 2M, but rise to of order 5% at the ends of the extended range 
M/5 < ji < 5M. The latter are fairly extreme scale choices, however. We believe that the 
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Figure 3: The CMS rapidity distribution of an on-shell Z boson at the LHC. The LO, NLO, and 
NNLO results have been included. The bands indicate the variation of the renormalization and 
factorization scales in the range M2/2 < /i < IMz- 



range used in the rest of the paper, [ip = fJ-R = n and M/2 < fj. < 2M, provides a good 
guide to the perturbative uncertainty remaining from the terms beyond NNLO. 

In Fig. 5 we present the rapidity distribution for on-shell Z production at Run II of 
the Tevatron. The scale variation is unnaturally small at LO; it is 3% at central rapidities, 
and varies from 0.1% to 5% from y = 1 to y = 2. This occurs because the direction of 
the scale variation reverses within the range of /i considered, i.e., daLo/dfi = for a value 
of fj, which satisifes Mz/2 < fj, < 2Mz. This value of /x depends upon rapidity, leading to 
scale dependences which vary strongly with Y. The scale variation exhibits a more proper 
behavior at NLO, starting at 3% at central rapidities and increasing to 5-6% at y = 2.5. 
At NNLO the scale dependence is drastically reduced, as at the LHC, and remains below 
1% for all relevant rapidity values. The magnitude of the higher-order corrections is slightly 
larger at the Tevatron than at the LHC. The NLO prediction is higher than the LO result 
by nearly 45% at central rapidities; this shift decreases to 30% at y = 1.5 and to 15% at 
y = 2.5. The NNLO corrections further increase the NLO prediction by 3-5% over the 
rapidity range y < 2. 

This remarkable stability of the rapidity distribution with respect to scale variation 
cannot be attributed to the smallness of the NNLO QCD corrections to the partonic cross 
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Figure 4: More general variations of the renormalization and factorization scales, for production 
of an on-shell Z boson at the LHC, at central rapidity Y = ^. For each order in perturbation 
theory (LO, NLO, NNLO), three curves are shown. The solid curves depict common variation of 
the renormalization and factorization scales, [ip ~ M_r as used in the rest of the paper, but 

extending the range of variation to A//5 < /i < hM. The dashed curves represent variation of the 
factorization scale alone, holding the renormalization scale fixed at M . The dotted curves result 
from varying the renormalization scale instead, holding the factorization scale fixed at M. 



sections. These corrections are the da'^^ jdY terms defined in Eq. (4.1) (after renormal- 
ization and mass factorization), convoluted with the MRST PDFs and with all partonic 
channels included. We vary the scale in these terms, and normalize this variation to the 
NLO cross section. We find that the NNLO corrections contribute a scale dependence 
of ~ 5% at central rapidities. When we form the complete NNLO cross section, which 
requires adding these corrections to the convolution of the da^^^ /dY and da^-^^ /dY terms 
of Eq. (4.1) with NNLO PDFs, the width of this band is decreased to less than 1%. This 
demonstrates a remarkable interplay between NNLO calculations and parton distribution 
functions. 

The small size of the NNLO corrections is partly due to large cancellations between 
the various partonic channels. To illustrate this, we present in Fig. 6 the fractional contri- 
butions of the various NNLO partonic corrections to the entire NNLO cross section, at Run 
I of the Tevatron. We include the qg and qiqj channels (the latter includes qq and qq inital 
states); the gg subprocess is numerically unimportant in this process. The magnitude of 
each order partonic correction, 6aij, can be 7-8% of the complete NNLO cross section. 
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Figure 5: The CMS rapidity distribution of an on-shell Z boson at Run II of the Tevatron. 
The LO, NLO, and NNLO results have been included. The bands indicate the variation of the 
renormalization and factorization scales in the range Mz/2 < n < 2Mz- 



cnnlO) at central rapidities, and can reach 10% of the entire result at larger rapidities. 
They cancel significantly, however, and their sum is only ~ 3% of the NNLO result. This 
cancellation is even larger at LHC energies; in fact, the qiqj and qg channels cancel to 
such an extent that the gg subprocess becomes an important contribution to the NNLO 
corrections. This split into partonic components is admittedly not entirely physical, as 
they are linked by initial-state collinear singularities. However, this degree of cancellation 
should be rather sensitive to the PDF set chosen. A different choice of PDFs may lead to 
changes in the cross section that are larger than that found by varying the renormalization 
and factorization scales. 

To investigate how the choice of PDFs affects the NNLO cross section, we first vary the 
MRST mode. The choices corresponding to the fast and slow DGLAP evolutions produce 
negligible shifts in our result, much less than 1% for all rapidities studied, and smaller than 
the residual scale dependence. (Similar results have been observed at the level of the total 
cross section [13].) However, the choice of MRST mode 4, which provides a better fit to 
the Tevatron high-E'r jet data, shifts the NNLO Z production cross section significantly. 
We present in Fig. 7 the rapidity distributions for LHC Z production using these two PDF 
choices. Both the NLO and NNLO results have been displayed; the scale variations are also 
included. The two mode choices are indistinguishable at NLO, due to the large residual 
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Figure 6: The fractional contribution of the various NNLO partonic channels to the entire NNLO 
cross section for Z production at Run 1 of the Tevatron for /i = Mz- Qiqj denotes all quark-quark 
and quark-antiquark channels, while qg indicates the quark-gluon and antiquark-gluon subprocesses. 
The gg channel is numerically small, and would be consistent with zero on this plot. 

scale dependence. At NNLO they become quite distinct, and the « 1% discrepancy is 
potentially visible given projected LHC errors. We note that the difference between the 
two PDF sets does not just produce a shift in the overall normalization. The mode 4 set 
slightly increases the number of quarks at x ~ 0.03, and decreases the number of gluons 
more substantially in this x range (to compensate for an even larger increase in g{x) at 
very large x). The qg channel has a negative partonic cross section; thus, paradoxically, 
decreasing g(x) increases the gluonic contribution to the cross section. The quark and 
gluon distribution shifts, plus a 2% increase in as{Mz), work in concert to increase the 
mode 4 predictions, relative to mode 1, for Z production at the LHC at central rapidities, 
and particularly in the range 1 < Y < 2. 

Another set of PDFs extracted with NNLO precision has been presented by Alekhin [14]. 
Only deep inelastic scattering data is used in this extraction; the NNLO QCD corrections 
can therefore be consistently included. The MRST global fits utilize processes for which 
these corrections are not known. This introduces an additional source of theoretical un- 
certainty into these parameterizations which is difficult to quantify. We present in Fig. 8 
a comparison between the MRST and Alekhin PDF sets for resonant Z production at the 
LHC. We have included the NNLO scale dependences for the Alekhin set and for the MRST 
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Figure 7: The rapidity distributions for Z production at the LHC for the MRST PDF sets mode 
1 and mode 4. The bands indicate the NNLO scale dependences, the solid lines denote the NLO 
mode 1 scale depedence, and the dashed lines indicate the NLO mode 4 scale variation. The upper 
lines correspond to the scale choice \x = 2M in the NLO cross sections, while the lower lines indicate 
^ = M/2. 

mode 1 and mode 4 sets; the NLO scale dependences for the MRST mode 1 and Alekhin 
parameterizations are also displayed. The large scale dependences again render all three 
choices indistinguishable at NLO. However, significant discrepancies appear at NNLO. The 
difference between the mode 1 and Alekhin sets is 2% at central rapidities; this increases to 
4.5% at y = 2 and to 8.5% at 1" = 3. The discrepancies in both normalization and shape 
will be clearly resolvable at the LHC. Although the MRST mode 4 choice is closer to both 
the shape and normalization of the Alekhin set, the differences still range from 1-8.5% as 
the rapidity is increased; this will again be observable at the LHC. Electroweak gauge bo- 
son production becomes a powerful discriminator between different PDF parameterizations 
when the NNLO QCD corrections are included. 

The di-lepton rapidity distribution for (Z, 7*) production has been measured by CDF 
at Run I of the Tevatron, in a mass window around Mz, 66 < M < 116 GeV [8]. To 
compare with these data, we numerically integrate over M as well as z and y in Eq. (6.3). 
The result is shown in Fig. 9. (The result of doing this M integral in a narrow-resonance 
approximation, taking into account the finite-mass endpoints, but neglecting photon ex- 
change, is about 2% lower.) The result with the Alekhin PDF set is about 4-5% above 
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Figure 8: The rapidity distributions for Z production at the LHC for the MRST PDF sets mode 
1 and mode 4, and for the Alekhin PDF set. The bands indicate the NNLO scale dependences; 
nil denoted the MRST mode 1 set, while ni4 indicates the MRST mode 4 set. The dashed lines 
denote the NLO scale dependence for the mode 1 set, and the dot-dashed lines denote the NLO 
scale dependence for the Alekhin set. The upper lines correspond to the scale choice /i = 2M in 
the NLO cross sections, while the lower lines indicate /i = Af/2. 



the MRST result. Naively, the Alekhin set gives a better fit to the data. However, most 
of the Alekhin/MRST difference here is in the overall normalization, and there is a 3.9% 
overall normalization uncertainty in the data (not shown in the error bars) due to the "pp 
luminosity uncertainty. Also, electroweak corrections have not yet been included. Hence 
the two PDF sets probably cannot be distinguished by this Run I data. Instead, it is clear 
from the figure that, for a given PDF set, the di-lepton rapidity distribution around the 
Z mass may be used to 'monitor' the luminosity at Run H, for which the statistical errors 
will be significantly smaller than those shown. 

We now examine the resonant production of W bosons at Run H of the Tevatron. We 
present in Fig. 10 the rapidity distribution for production; the distribution for the W~ 
can be obtained by substituting Y — > —Y. Both the scale variations and the magnitudes 
of the higher corrections are similar to those found previously for Z production at the 
Tevatron. The scale dependence at LO is again unnaturally small, ranging from 3-5%, 
because daLo/dn = for values of within the parameter space studied. At NLO the scale 
variations are between 2% and 3.5%; they decrease to 0.3-0.7% at NNLO, depending 
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Figure 9: The di-lepton rapidity distribution for (Z, 7*) production at Run I of tlie Tevatron, 
compared witli data from CDF [8]. Tlie LO and NLO curves are for tlie MRST PDF set. Tlie 
thin NNLO bands are for the MRST (lower) and Alekhin (upper) parameterizations. The bands 
correspond to varying M/2 < < 2M . 



upon the rapidity chosen. The magnitude of the NLO corrections is large, varying from 
45% at central rapidities to ~ 25% at larger rapidities. The NNLO corrections are also 
appreciable; they range from 2.5% at 1" = to 4% at \Y\ ~ 2. 

Another observable frequently studied at hadron colliders is the W charge asymmetry, 
defined as 

^^^_ da{W+)/dY-da{W')/dY 
^^^^^ " da{W+)/dY + da{W-)/dy ^^'^^ 

A simple calculation in the LO approximation reveals that this quantity is sensitive to the 
X dependence of u(x)/d{x), the ratio of up and down quark distributions in the proton. 
Although in a realistic experiment only the pseudorapidity of the charged lepton coming 
from the W decay can be measured, much of the sensitivity to the PDFs remains. Since 
Aw is a ratio of cross sections, it might be expected that it is rather insensitive to QCD 
corrections. This is indeed the case. At the Tevatron, a pp collider, with the assumption 
of CP invariance, the charge asymmetry is an odd function of Y, since it may be written 
as 

A (Y\ da{W+)/dY-da{W+)/dY\Y^-Y _ , , 

= da{W+)/dY + da{W+)/dY\y^^y " -^^(-^)- (^.8) 
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Figure 10: The CMS rapidity distribution of an on-shell W'^ boson at Run II of the Tevatron. 
Shown are the LO, NLO, and NNLO results for the MRST PDF sets. The bands indicate the 
variation of the renormahzation and factorization scales in the range Mw/'2 < < 2Mw- 



The asymmetry is positive for positive Y, corresponding to the boson moving in the 
same direction as the incident proton, because u{x) is larger than d{x) at large x. In 
Fig. 11, we present the LO, NLO, and NNLO predictions for the charge asymmetry at Run 
II of the Tevatron, together with their scale dependences. The NLO corrections increase 
the Born level result by 2-4%. The NNLO corrections to the NLO result range from —2% 
at central W rapidities to +1% at large rapidities. 

The scale variations of Aw{Y) are small; to study them, we present in Fig. 12 the 
scale-dependence bandwidths, defined as 

AwjY, ^ = 2Mw) - Aw{Y, = Mw/2) 

The scale variation is already below 5% for all rapidities at LO, and is below 1% at NLO. 
The NNLO prediction is absolutely stable against scale variation, indicating that this 
observable is potentially a very strong constraint on quark distribution functions. We note 
that the scale choice fi = 2Mz in the LO asymmetry yields an approximation to the NNLO 
result which is accurate to 1-2% for essentially all rapidities. 



The NNLO predictions for the rapidity distributions for on-shell W boson produc- 
tion at the LHC are shown in Fig. 13. The distributions are symmetric in Y; only the 
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Figure 11: The W charge asymmetry at Run II of the Tevatron. Included are the LO, NLO, and 
NNLO results. The bands indicate the variation of the renormalization and factorization scales in 
the range Mw/'2 < H < 2M\y- As the charge asymmetry is rather insensitive to QCD corrections, 
the three bands are almost completely degenerate. 

positive half of the rapidity range is shown for , and the negative half for W~. The 
charge asymmetry is positive for all rapidities, but is particularly striking around 1" = 3. 
The behavior of the perturbation series is very similar to that discussed previously for Z 
production at the LHC. Again the NNLO scale- variation bandwidths are extremely narrow 
for central rapidities, ranging from « 0.6% for 1" < 2, to 1.5% at 1" = 3, to 3% at y = 4. 

In addition to the study of resonant production of electroweak gauge bosons, both 
the Tevatron and the LHC use high-invariant-mass Drell-Yan production of lepton pairs 
to search for new gauge bosons and lepton-quark contact interactions. Although these are 
primarily inclusive searches, rapidity cuts are required because of experimental constraints. 
We therefore examine the NNLO QCD corrections to off-shell (Z, 7*) production at large 
invariant masses. We present below the rapidity distribution for M = 250 GeV (Z, 7*) 
production at the LHC in Fig. 14, and for M = 200 GeV at Run II of the Tevatron in 
Fig. 15. The scale dependences are significantly smaller for M = 250 GeV than for resonant 
Z production at the LHC. The LO scale variation is 12% at central rapidities and 4% at 
Y = 3. Both the NLO and NNLO scale variations are much less than 1% for all values 
of rapidity. The magnitude of the higher-order corrections is much larger, however. The 
NLO result increases the LO prediction by nearly 35% at central rapidities; this correction 
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decreases to 10% at larger Y values. This discrepancy between the sizes of the scale 
variations and of the NLO shifts sends a somewhat mixed message regarding the importance 
of the NNLO corrections. We find that they are small, decreasing the NLO result by less 
than 0.5% for Y < 1.5, and increasing it by less than 1% for 1.5 <Y< 2.8. The small scale 
dependence of the NNLO cross section and the stability of the NLO prediction indicate a 
complete stabilization of the perturbative result for M = 250 GeV at the LHC. 

The results for M = 200 GeV (Z, 7*) production at Run II of the Tevatron exhibit 
both larger scale dependences and more important higher-order corrections. The LO scale 
variations are similar to those found at the LHC, ranging from 7% at y = to ~ 15% at 
larger rapidity values. In contrast to the LHC case, the NLO scale dependences remain 
fairly large, varying from 5% at central rapidities to 14% at 1" = 2. At NNLO, the scale 
variations are between 1.5% and 4%, again increasing for larger rapidities. The magnitude 
of the NLO corrections is over 40% at central rapidities, and « 30% at larger Y values. 
The NNLO corrections further increase the NNLO result by 5-6% throughout the entire 
rapidity range. 



Finally, we study the accuracy of various approximations to the complete NNLO cor- 
rection to the rapidity distribution. There are three distinct types of terms which appear 
in the result: 
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Figure 13: The CMS rapidity distributions for production of an on-shell W~ boson (left) and 
on-shell W+ boson (right) at the LHC, at LO, NLO, and NNLO, for the MRST PDF sets. Each 
distribution is symmetric in Y; we only show half the rapidity range in each case. The bands indicate 
the common variation of the renormalization and factorization scales in the range Mw/'2 < fJ- < 
2Mw 

• soft (sz)'- terms which contain either a delta function or a plus distribution in 1 — 
z. These terms arise from production of the vector boson V close to the partonic 
threshold, and can be obtained by considering only soft partonic emissions from the 
qq ^ V subprocess. 

• coUinear (cy): terms containing delta functions or plus distributions in either y or 
1 — y, but not in 1 — z. These terms result from the emission of radiation collinear 
to one of the initial partons. 

• hard {h): terms which have no delta functions or plus distributions. These terms 
arise from generic scattering events with the emission of hard additional partons in 
the final state. 

There is some potential ambiguity in this separation, due to the presence of Jacobian 
factors in the integration. We perform the separation in terms of the functions Fij{z,y) 
appearing in Eq. (6.3); i.e., including all Jacobian factors resulting from the transformation 
the variables {z,y). The Sz terms can be obtained by using the soft gluon approximation, 
and it is possible to imagine obtaining the Cy contributions from a simplified calculation 
in which the collinear emission of V is factorized from a hard scattering piece. The hard 
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Figure 14: The rapidity distribution for (Z, 7*) production at the LHC for an invariant mass 
M = 250 GeV. The LO, NLO, and NNLO results have been included. The bands indicate the 
residual scale dependences. 

emissions, however, require a full NNLO computation. Intuitively, we expect the Sz terms, 
which are the simplest to obtain, to dominate for large invariant masses, i.e., as the z — > 1 
threshold is approached. We wish to examine whether this contribution, or perhaps the 
Sz and Cy terms together, can furnish a reasonable approximation in phenomenologically 
interesting regions of parameter space. 

We present in Figs. 16 and 17 the NNLO corrections to the rapidity distributions 
for (.2', 7*) production at the LHC, split into its soft, collinear and hard components, for 
the invariant masses M = Mz and M = 2 TeV. The NNLO corrections are the da^'^^/dY 
terms defined in Eq. (4.1), convoluted with the MRST PDFs and with all partonic channels 
included. We present separately the following pieces: the Sz term, the Cy term, the h term, 
and the sum of the h and Cy pieces, which would integrate to the "hard" (non-soft) part 
of the total cross section. These terms are normalized to the complete NNLO correction. 
At M = Mz, all components are important. We note that there are large cancellations 
between the Sz term and the remaining pieces. Neither the Sz piece nor the sum of the 
Sz and Cy terms furnishes a good approximation to the complete result. Generic hard 
emissions are important; this result is expected, since there is a large amount of phase 
space available. At M = 2 TeV, the magnitude of the Sz term becomes larger compared to 
the hard and Cy terms, as expected. However, it still does not furnish a good approximation 
to the entire result for all rapidities; the fact that it does so for central rapidities arises from 
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Figure 15: The rapidity distribution for (Z, 7*) production at Run II of the Tevatron for an 
invariant mass M = 200 GeV. The LO, NLO, and NNLO results have been included. The bands 
indicate the residual scale dependences. 



an accidental cancellation between the hard and Cy pieces. We observe similar behavior for 
Tevatron kinematics. We note that at higher invariant masses, the magnitude of the hard 
term decreases quickly. The Cy term also decreases, but less rapidly. The Sz term does not 
dominate until very large invariant masses are reached. 



7. Conclusions 

We have presented a calculation of the rapidity distributions for electroweak gauge boson 
production at hadron colliders through NNLO in QCD. This is the first complete NNLO 
computation of a differential quantity needed for high-energy hadron collider physics. We 
have discussed in detail a powerful new technique for calculating differential distributions. 
This method is completely automated, produces fully analytic results, and treats the vari- 
ous components of a NNLO calculation in a unified manner. Our results will assist in the 
extraction of parton distribution functions, parton-parton luminosities, electroweak gauge 
boson information, and other quantities of interest with the accuracy needed for Tevatron 
and LHC physics. 
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Figure 16: The components of the NNLO corrections to the rapidity distribution for (Z, 7*) 
production at the LHC for M = Mz- The pieces included are the hard part /i, Sz, Cj,, and the sum 
of the h and Cy pieces. The complete NNLO correction h + Cy-\- Sz is normalized to unity. We have 
set /I = Af . 



We have found that the residual scale dependences for resonant W and Z produc- 
tion at both the Tevatron and the LHC are below 1% when the NNLO corrections are 
included; the rapidity distributions are completely stable against higher-order QCD cor- 
rections. Only higher-order electroweak corrections and mixed QCD-electroweak effects 
remain to be included [11]. These distributions are therefore ideal observables to use to 
discriminate between different parton distribution function parameterizations. We have 
studied several different NNLO extractions of parton distribution functions obtained by 
the MRST group, as well as an NNLO extraction provided by Alekhin. Varying the evolu- 
tion rate of the approximate NNLO DGLAP kernels in the MRST parameterization yields 
negligible shifts in our results. However, an MRST PDF set designed to provide a better 
fit to the Tevatron high-i?T jet cross section produces a difference of about 1% in rapidity 
distributions at the LHC. This difference may be observable, given expected experimental 
errors. 

The deviations induced by instead using Alekhin's PDF extraction are more striking. 
Both the normalization and shape of the rapidity distributions obtained with Alekhin's 
parameterization differ from those found with the MRST sets; the differences range from 
2-8.5% as the rapidity is varied. These differences should be easily resolvable at the LHC, 
given the expected errors. The MRST parameterizations are derived from global fits to a 
variety of data, including data from processes for which the NNLO QCD corrections are 
unknown. We note that the magnitude of the discrepancies between the Alekhin and MRST 
PDF sets is consistent with the typical size of NNLO QCD corrections. It is conceivable that 
the inclusion of these corrections into the MRST fit might lessen the observed differences. 
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Figure 17: The components of the NNLO corrections to the rapidity distribution for (Z, 7*) 
production at the LHC for M = 2 TeV. The pieces included are the hard part /i, s^-, Cy, and the 
sum of the h and Cy pieces. The complete NNLO correction h + Cy + Sz is normalized to unity. We 
have set fi = M. 



In fact, the NNLO Alekhin PDF set includes a full error matrix. (Similar uncertainty 
estimates are available for the MRST set at NLO.) This matrix permits the construction 
of PDF uncertainty bands for the vector boson rapidity distributions, whereas here we just 
employed the central PDF values. We defer such a study to future work. 

The magnitude of the NNLO corrections to resonant gauge boson production ranges 
from 1-2% at the LHC to 3-4% at the Tevatron; the corrections for higher invariant mass 
gauge bosons can reach 5-6% at the Tevatron. These contributions must be included 
to yield a theoretical calculation accurate to ~ 1%, the projected experimental precision 
at the LHC. However, the NNLO corrections do not vary strongly with rapidity. The 
NLO rapidity distribution appears to describe the kinematics quite well. Reweighting the 
NLO distributions by the inclusive K-factor K^'^'> = o"nnlo/c"nlo yields an approximation 
accurate to < 1% for all relevant rapidities. The analogous reweighting of the LO results, by 
^NNLO _ (7|v^f^Lo /clo ) does not furnish a good approximation to the complete result. The 
excellent accuracy of the NLO reweighting technique for the rapidity distribution suggests 
that one apply the factor K^"^^ to output from a hadron-level Monte Carlo program which 
incorporates the NLO vector-boson production matrix elements, such as MC@NLO 2.2 [54]. 
This simple procedure should give a good picture of the structure of the hadronic events 
accompanying the vector bosons, and is likely to approach NNLO precision for sufficiently 
inclusive observables. 

We have also studied the accuracy of approximating the NNLO corrections by partial 
results. We have found that including only virtual and soft gluon corrections, labeled as 
Sz in the text, does not yield a good approximation for resonant gauge boson production. 
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Only at very large invariant masses do these terms dominate. We estimate that average 
values of Bjorken x > 0.3 — 0.4 must be reached before the Sz component accounts for 
80% of the complete NNLO correction for all relevant rapidities. We also note that the 
Sz terms do not accurately predict the shape of the NNLO correction, as is apparent from 
Figs. 16 and 17. 

Finally, we note that with our result for the rapidity distribution, it is possible to 
obtain almost full control over the kinematics of the electroweak boson, as produced in 
fixed-order perturbation theory. This is because the NLO QCD corrections to the double 
differential distribution (Pa/{dYdp±) for electroweak boson production are known [55]. It 
was assumed in Ref. [55] that p± ^ 0. It is therefore not possible to perform the integration 
over p_i to get da/dY using the results of Ref. [55] alone. However, the NNLO calculation 
of the rapidity distribution presented here gives an unambiguous answer for the integral 
over p_\_ at fixed values of rapidity, and can therefore be used as a normalization condition. 
We write 

e(pT-pi), (7.1) 

where d'^a/(dYdp±) is the distribution computed in Ref. [55]. Integrating d'^a„io(\/{dYdp±) 
over p± gives the correct result for the rapidity distribution; however, the "zero p±" bin 
extends from pj_ = to p± = ]9^*. Apart from this drawback, Eq. (7.1) provides a simple 
way to describe the electroweak boson kinematics at NNLO in QCD. 

Our results are an important theoretical input for physics at both the Tevatron and 
the LHC. We believe the method we have introduced to obtain these results can be used to 
calculate other phenomenologically interesting observables. We anticipate its application 
in many other areas of collider physics. 
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